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Abstract The propagation of solar waves through the sunspot of AR9787 is ob- 
served using temporal cross-correlations of SOHO/MDI Dopplergrams. We then 
use three-dimensional MHD numerical simulations to compute the propagation 
of wave packets through self-similar magneto-hydrostatic sunspot models. The 
simulations are set up in such a way as to allow a comparison with observed 
cross-covariances (except in the immediate vicinity of the sunspot). We find 
that the simulation and the /-mode observations are in good agreement when 
the model sunspot has a peak field strength of 3 kG at the photosphere, less so 
for lower field strengths. Constraining the sunspot model with helioseismology 
is only possible because the direct effect of the magnetic field on the waves has 
been fully taken into account. Our work shows that the full- waveform modeling 
of sunspots is feasible. 
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1. Introduction 

The subsurface magnetic structure of sunspots is poorly known. The theoretical 
picture is that sunspots are either monolithic or not (Parker, 1979), or change 
from being monolithic to "disconnected" over the course of the first few days 
of their lives (Schussler and Rempel, 2005). The observational picture is limited 
because the subsurface structure is inherently difficult to infer. The most promis- 
ing possibility by far is local helioseismology. Local helioseismology includes 
several techniques of data analysis, such as Fouricr-Hankel analysis (e.g. Braun, 
1995), time-distance analysis (e.g. Duvall et ai, 1993) and helioseismic holog- 
raphy (e.g. Lindsey and Braun, 1997). For example the time-distance approach 
(Duvall et al, 1993) has been applied to determine wave-speed variations and 
flows associated with sunspots (e.g. Kosovichev, Duvall, and Scherrer, 2000; 
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Zhao, Kosovichev, and Duvall, 2001; Couvidat, Birch, and Kosovichev, 2006). 
These inversions did not take into account the direct effects of the magnetic 
field. We do not mean to give the impression that such effects have not been 
investigated, but rather that they are only beginning to be incorporated into 
helioseismic inversions. Fourier-Hankel decompositions of the acoustic wave field 
(p modes) near sunspots by Braun, Duvall, and LaBrontc (1987) showed that 
incoming p modes arc phase-shifted and "absorbed" by sunspots. This triggered 
many studies of the effects of the magnetic field on solar waves. For example, 
Spruit (1991) suggested that sunspot magnetic fields are responsible for the 
observed acoustic wave "absorption" by partially converting incoming p modes 
into slow magnetoacoustic waves. This idea was followed up in detail by Spruit 
and Bogdan (1992), Cally and Bogdan (1993), Cally, Bogdan, and Zweibel 
(1994), Hindman, Zweibel, and Cally (1996), Cally and Bogdan (1997), Bogdan 
and Cally (1997), and Rosenthal and Julicn (2000). An important result was 
the realization that non-uniform and non-vertical magnetic fields are required 
to explain the observations. In particular, Cally (2000) reported on numerical 
calculations showing that inclined magnetic fields arc able to achieve levels of 
absorption compatible with the observations. This aspect of the problem was 
followed up by Cally, Crouch, and Braun (2003) and Crouch et al. (2005), which 
placed constraints on the strength of the sunspot's magnetic field. More recently 
attention has been placed on upward propagating magnetoacoustic waves and 
their possible observational signatures (Schunkcr and Cally, 2006; Khomenko 
and Collados, 2006; Cally, 2007; Cally and Goosscns, 2008). The use of direct 
numerical simulations in helioseismology, with or without magnetic fields, has 
also undergone rapid development. The aims of these simulations include val- 
idating helioseismic techniques (e.g. Zhao et al., 2007) and, more importantly, 
increasing our understanding of what the seismic observations reveal about the 
solar interior. Here we restrict our focus to simulations of linearized wave prop- 
agation through an inhomogeneous solar atmosphere. Examples of these linear 
wave propagation studies, most of which do not include the direct effects of 
magnetic fields, include Birch et al. (2001), Khomenko and Collados (2006), 
Hanasoge, Duvall, and Couvidat (2007), Hanasoge et al. (2006), and Hanasoge 
and Duvall (2007). Some results of these studies clearly have a large bearing on 
the correct interpretation of the helioseismic signatures of sunspots. 

In this paper we first design a technique to image the interaction of waves with 
sunspots, using appropriate averaging of the cross-correlations of the random 
seismic wave field (SOHO/MDI data). We proceed further by using numerical 
simulations to do full-waveform forward modeling of the passage of solar waves 
through a sunspot. This is done with the SliM code (Cameron, Gizon, and 
Daiffallah, 2007), specifically developed for this purpose. The ultimate goal is 
to understand the observed cross-correlations in terms of the properties of our 
parametric sunspot models. 

This paper is organized as follows. The SOHO/MDI observations of sunspot 
AR9787 are presented in Section 2 and their helioseismic analysis in Section 3. 
We describe our numerical simulation code in Section 4. We compare the obser- 
vations (cross-correlations) against the vertical velocity on a horizontal cut taken 
from the simulation at the height of the quiet-Sun photosphere. This, amongst 
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other reasons, means the comparison is only meaningful outside the sunspot. 
Except in the immediate vicinity of the sunspot, the comparison between the 
observations and the simulations is very encouraging, as shown in Section 5. We 
place a helioseismic constraint on the sunspot magnetic field in Section 6. Our 
work strongly suggests that we will be able to use observations and simulations 
in combination to constrain the subsurface structure of sunspots by taking direct 
and indirect magnetic effects into account. 

2. SOHO/MDI Observations of Active Region 9787 

The helioseismic analysis of a sunspot is more or less difficult depending on which 
sunspot is being studied. We have been searching the SOHO/MDI database for 
the ideal sunspot, i.e. a sunspot which is isolated (does not belong to a complex 
sunspot group), has a simple geometry (circular shape), and evolves slowly in 
time as it moves across the solar disk. In order to find such a "theorist's sunspot" , 
we searched the MDI Dynamics Program which combines the advantage of a 
good spatial resolution (0.12° at disk center) with a complete view of the solar 
disk. The full-disk Dopplergrams are available each minute during two to three 
months each year since 1996. The sunspot we have selected is that of Active 
Region 9787, continuously observed by MDI during nine days: 20-28 January 
2002. The Dopplergrams were remapped using Postel's azimuthal equidistant 
projection almost - but not exactly - centered on the sunspot, using a tracking 
angular velocity of —0.1102° day -1 (in the Carrington frame). In addition to 
the Dopplergrams, we also used the line-of-sight magnetograms (each minute) 
and all the intensity images (one per six hours). The daily averages of these 
three quantities arc shown in Figure 1. Apart from some plage, there is no 
other active region in the vicinity of the sunspot during the entire observation 
sequence. The sunspot of AR9787 is large and quite stable over the nine days 
of the observations, even though it starts decaying from 27 January onward. 

Using the MDI intensity images, we measured the center of the sunspot's 
umbra. As shown in Figure 2, the sunspot has a significant amount of proper 
motion. Thus we chose to split the time series into six-hour subsets and to 
analyze each subset separately. The images belonging to each particular subset 
were remapped into a new Postcl map, with the center of the projection cor- 
responding exactly to the center of the sunspot's umbra in the middle of the 
six-hour time interval. These 36 six-hour time scries of Dopplergrams, centered 
on the sunspot's position, are the basic data that we used for the helioseismic 
analysis, which we discuss in the next section. 

The average intensity image of the sunspot, corrected for the proper motion 
of the sunspot, is shown in Figure 3. From this average image, which is still 
sharp, we determine the average umbral and penumbral radii to be 9 Mm and 
20 Mm respectively. 

3. Observed /-mode Cross-Covariance Function 

We study sunspot AR9787 using cross-covariances of the random wavefield, as 
is done in time-distance helioseismology (Duvall et al, 1993). The temporal 
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Figure 1. SOHO/MDI observations of the sunspot of Active Region 9787 during the period 
20 — 28 January 2002. Shown are the daily averages of the line-of-sight Doppler velocity (top 
row), the intensity (middle row), and the line-of-sight magnetic field (bottom row). The color 
bars are in units of kms -1 , relative intensity, and kG, from top to bottom. 
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Figure 2. Coordinates of the center of the sunspot's umbra as a function of time. Left: 
Carrington longitude of the umbra (dots) and Carrington longitude at the center of the Postel 
projection (line segments). Right: Latitude of the umbra (dots) and latitude at the center of 
the Postel projection (horizontal line). 



cross-covariance between two points on the solar surface provides information 
about the Green's function between these two points. This interpretation has 
recently been shown to be correct in the case of homogeneously distributed 
random sources in an arbitrarily complex medium (Colin de Verdiere, 2006, 
Gouedard et al, 2007). 

In this paper, for the sake of simplicity, we wish to observe the propagation 
of /modes through the sunspot. Thus we filter the Dopplergrams in 3D Fourier 
space to only keep the / modes. Let us denote by 4>(r,t) the filtered Doppler 
velocity, where r is a position vector and t is time. Rather than studying the 
cross-covariance of <f> between two spatial points, we consider the cross-covariance 
between <p averaged over a great circle 7, denoted by ^(7, t), and <j) measured at 
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Figure 3. SOHO/MDI intensity image of AR9787 averaged over nine days, after correcting 
for the proper motion of the sunspot. The intensity is measured relatively to the quiet-Sun 
value. 



any other point r: 



where the effective integration time is T is nine days, or 36 times six hours. Since 
averaging <j) over 7 is equivalent to selecting only the horizontal wavevectors that 
are perpendicular to 7, the above cross-covariance function gives us information 
about plane wave packets that propagate away from 7. In simple words, the 
cross-covariance at time lag t tells us about the position of a wave packet, a time 
t after it has left 7. In the rest of the paper we fix the distance between 7 and 
the center of the sunspot at A = 40 Mm. Since A>A, where A ~ 5 Mm is the 
dominant wavelength, the sunspot is in the far field of 7. 

For any particular choice of orientation of 7, the computed cross-covariance is 
very noisy. Thus the need for some spatial averaging. Let us pick a reference great 
circle coincident with the meridian at a distance 40 Mm eastward of the center of 
the sunspot. Because sunspot AR9787 is almost rotationally invariant around its 
center, we can compute many equivalent cross-covariance functions correspond- 
ing to many different directions of the incoming wavepacket, derotate these about 
the center of the sunspot so that they match the reference cross-covariance, and 
average them to reduce the noise. We have performed this averaging over all the 
directions of the incoming wavevectors, with a fine sampling of 1°. As seen in 
Figure 4, this enables us to reach a very good signal-to-noise level. 
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Figure 4. The upper panel shows the observed /-mode MDI cross-covariance function at zero 
time lag. The sunspot (black circle) is at the center of the Postel map, a distance A = 40 Mm 
from the meridian 7 (white). The lower panel shows the same cross-covariance at time lag 
t = 130 minutes. In order to increase the signal-to-noise ratio, the cross-covariance was averaged 
over all wave directions. It is easy to see, by eye, that the waves are affected by their passage 
through the sunspot. The white rectangles show the size of the simulation box. 
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4. Three-dimensional MHD Simulation of Wave Propagation 
through a Sunspot 

4.1. The Equations 

We use the ideal MHD equations linearized about an arbitrary, inhomogeneous, 
magnetized atmosphere. We assume a local Cartesian geometry defined by hori- 
zontal coordinates x and y, and the vertical coordinate z increases upwards. The 
level z = is assumed to correspond to the photosphere, as defined in model S 
(Christensen-Dalsgaard et al, 1996). Under the assumptions that gravitational 
acceleration (g) is constant and that there is no background steady flow, the 
equation governing the wave- induced displacement vector £ is (e.g. Cameron, 
Gizon, and Daiffallah, 2007) 

PodU = F', (2) 

where 

F' = - VP' + p'gz + -!- (J' x B + J x B') (3) 

47T 

is the linearized force acting on a fluid element (first order in £). In this paper 
the subscript variables represent the steady inhomogeneous background atmo- 
sphere and the primed quantities represent the wave-induced perturbations. In 
the above equation, the density, pressure, magnetic field, and electric current are 
denoted by the symbols p, P, B and J respectively. The system is closed by the 
following relations that define F' in terms of £: 



P = -v ■(/*,€), (4) 

p' = coV + rvpo)-rvPo, (5) 

B' = V x (£ x B ), (6) 

J' = V x B', (7) 



where Co is the background sound speed. 

Whilst we have employed ideal MHD, the waves on the Sun are strongly atten- 
uated, presumably as a result of scattering off the time-dependent granulation. 
Empirically, a solar mode with horizontal wavevector k decays as e~ lkt where 
l/7fe is the e- folding lifetime at wavenumber k — \\k\\. Introducing v' as the wave 
velocity, for each horizontal Fourier mode we write 

Po(d t +jkW(k,z,t) = F'(k,z,t), (8) 
(9t+7fc)€(M,*) = v'(k,z,t). (9) 

In the above equations, any function f(k, z, t) refers to the spatial Fourier trans- 
form of f(x, y, z, t). For 7fe > the system is damped. The eigenstates, pairs of 
j arc independent of 7^, although the eigenvalues are naturally affected. 
This is the advantage of this approach. The attenuation that we have imple- 
mented acts only in the time domain and is suitable for our purpose, although 
in the Sun attenuation is more complicated. An alternate way of introducing 
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this phenomenological time attenuation would have been to perform the ideal 
calculation and then impose the decay in the time domain post-facto after the 
calculation has ended. 

In this paper we focus on the / modes (solar-surface gravity waves), for which 
the attenuation has been measured to be 

Ik = 7* (fc/fc*) 2 ' 2 , (10) 

where j*/tt = 100 /iHz and fc* = 902/_R© is a reference wavenumbcr (Duvall, 
Kosovichev, and Murawski, 1998; Gizon and Birch, 2002). 

4.2. The Code 

We use a modified version of the SLiM code which is, apart from the modifica- 
tions discussed below, described in Cameron, Gizon, and Daiffallah (2007). The 
code has been tested against analytic solutions - some of these tests are also 
described in detail in Cameron, Gizon, and Daiffallah (2007). 

The present code includes two absorbing layers at the top and the bottom of 
the box. In the top layer above the temperature minimum we heavily damp the 
waves and systematically reduce the effect of the Lorentz force. This layer only 
affects waves which have escaped through the photosphere. Likewise the bottom 
layer damps the waves that propagate downwards. The purpose of both layers 
is to minimize the effects of the boundaries which would otherwise artificially 
reflect the waves. 

An additional change has been to introduce the mode attenuation described 
in Section 4.1. Since we use a semi-spectral scheme the implementation was 
straightforward. 

The scheme uses finite differences in the vertical direction, with 558 uniformly 
spaced grid points sampled at Az — 25 km. In the horizontal direction we use a 
Fourier decomposition with 200 modes in the x direction and only 50 in the y. 
The spatial sampling is Ax = 0.725 Mm and Ay = 1.45 Mm. This seemingly low 
resolution is satisfactory because neither the initial wave-packet nor the sunspot 
has any significant power at short wavelengths. The size of the simulation box 
is 145 Mm long (x-coordinate), 72.5 Mm across (^-coordinate) and 14 Mm in 
depth (12.5 Mm below the photosphere). A typical run, such as the one from 
Section 5.3, takes approximately 14 days on a single CPU core. 

4.3. Stabilizing the Quiet-Sun Atmosphere 

Our aim is to model the propagation of waves through the solar atmosphere in 
such a way as to allow a direct comparison with the observations. The most direct 
way of proceeding would be to use an existing solar-like atmosphere such as that 
of model S (Christensen-Dalsgaard et al, 1996). This most direct approach is 
however unavailable when considering the full evolution of wavepackets using nu- 
merical codes because both the solar and model S atmospheres are convectively 
unstable. 

Wave propagation through unstable atmospheres is difficult to study numeri- 
cally because any numerical noise in the unstable modes grows exponentially and 
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eventually dominates the numerical solution. The problem is then to stabilize 
the atmosphere whilst leaving it as solar-like as possible, at least in terms of the 
nature of the waves propagating through it. 

The convective instability is easily understood by reference to Equation (5). 
We imagine a small vertical displacement of a blob of plasma, and assume it 
is evolving slowly enough that it is in pressure balance with its surroundings. 
This implies P' = and hence CQp' = £ • VP — Cq£ • Vp - For a vertically 
stratified atmosphere, this becomes p' = £, z (d z Po — c o^Po)/cq- The atmosphere 
is convectively unstable if an upward displacement (£ z > 0) corresponds to a 
region of lowered density (p' < 0) ; in such a case the fluid parcel is buoyant and 
accelerates upwards. The atmosphere is convectively stable when £ 2 and p 1 have 
the same sign, requiring (d z Po — Cp9 z po)/ c o > 0; equivalently 

d z P a > c 2 d z p . (11) 

We have the freedom to modify any combination of P , c , or p in order to 
satisfy Equation (11). It is somewhat natural to regard these three quantities as 
being related through an equation of state or through the constraint that the 
atmosphere be hydrostatic, however neither of these relationships is necessary 
in terms of the properties of the propagating waves. We choose Po, cq, and po 
so that the wave speed is solar-like. Changing the sound speed would obviously 
have a major impact on the propagation of sound waves, and hence we have 
chosen to keep Co unchanged. Varying po would have the effect of varying the 
distribution of the kinetic energy density of the different modes as a function 
of height, which would significantly change the sensitivity of wave packets to 
inhomogeneities; so we have also kept p fixed. Thus we choose P so that 

d z P Q = max{O.9cg0*po, d z P u }, (12) 

where P u was the pressure distribution of the unstable atmosphere. Notice that 
d z Po is negative so that the factor of 0.9 does indeed mean that Equation (11) is 
satisfied. The factor of 0.9 is possibly unnecessary, but does little harm: setting 
this factor to 1 should create stationary eigenmodes, setting it to 0.9 introduces 
internal gravity modes which propagate very slowly. 

4.4. An Example Oscillation Power Spectrum 

Since we have modified the atmosphere, we need to check whether it supports 
oscillations that are those of the Sun and model S in the absence of a sunspot. 
We chose initial conditions consisting of a "line source" of the form 

& = e -x 2 /2s 2 e -(z-zo) 2 /2s 2 att = 0, (13) 

where s = 0.7 Mm and z = —250 km is an arbitrary height below the photo- 
sphere. All other wave perturbations, £ x , £ y , and v' being zero at t = 0. This 
disturbance encompasses many different modes allowing several ridges of the 
dispersion diagram to emerge. Figure 5 directly compares these ridges with the 
eigenfrequencies of model S (A.C. Birch, private communication). The / modes 
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Figure 5. Power as a function of wavcnumber and frequency for the example simulation 
described in Section 4.4 (no sunspot). The dashed lines correspond to the eigenfrequencies of 
the model S atmosphere. 



and the acoustic modes pi, p2, and P3 lie where they are expected for kR Q > 300. 
It is also to be noted that there is a very low amplitude ridge (not visible on the 
plot) at frequencies below 1 mHz that corresponds to internal gravity modes, 
which are not present in the Sun. They are an artifact of having made the system 
convectively stable as explained above, but they can be safely ignored. 

The simulations and model S eigenfrequencies differ strongly at wavenumbers 
less than 3OO/i?0. This difference is due to the fact that our computational 
domain extends to only 12 Mm below the solar surface. We have not tested a 
deeper box because the /, p\, and p 2 modes look satisfactory in the range of 
wavelengths we are interested in (k > 300/i? Q ). 
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4.5. A parametric sunspot model 

In this section we will describe the sunspots that we have embedded in our 
atmosphere. We begin by noting that we embed the sunspot in the atmosphere 
before we stabilize it, in order to get the correct sound speed and density 
structures. We use P u to denote the pressure of the atmosphere before it has 
been stabilized as described above. The sunspot is modeled by an axisymmetric 
magnetic field B = B (r, z) , where r is the horizontal radial distance from 
the sunspot axis. Within the part of the atmosphere described by model S, 
the magnetic field is made hydrostatic in a standard way, by calculating the 
Lorentz force L = (J x B )/47r and noting that the horizontal force balance 
then requires a horizontal pressure gradient: d r P u = f • L. Since P u is unaffected 
by the spot at large distances, we can integrate from infinity towards the center 
of the spot to find Po- Having thus found P u we can find p by the constraint 
of vertical force balance. In this case the gravitational force needs to balance 
both the vertical component of the Lorentz force as well as that of the pressure 
gradient. In principle given p and Po the Saha-Boltzmann equations need to 
be solved to obtain the first adiabatic exponent (ri) and the sound speed (co). 
However for the purposes of this paper we have assumed that Ti is a function 
of z only and is unaffected by the sunspot. This assumption will be relaxed in 
future studies. 

The procedure thus outlined can always be applied, however in some cases the 
results will involve negative pressures or densities. Such solutions are of course 
unphysical and indicate that no hydrostatic solution exists for the given magnetic 
configuration and quiet-Sun pressure stratification. The problem typically arises 
in the very upper layers of the box when the magnetic pressure is large. For 
example a purely vertical flux tube with a magnetic pressure P m cannot be 
in hydrostatic balance in an atmosphere with external pressure P cxt < P m . In 
practice the Sun rapidly evolves to almost force-free field configurations above 
the solar surface. The role of the Lorentz force in structuring the atmosphere 
in the low-/? region is thus artificial as well as being responsible for the lack of 
equilibrium. It is also dynamically of minor importance to the waves since it is 
above both the acoustic cut-off and the layer where the acoustic and Alfven wave 
speeds become equal. We have therefore adopted the approach (A. Nordlund, 
in the context of realistic photosphcric magnetoconvection simulations of active 
regions, private communication) of scaling the Lorentz force in this region when 
constructing the background atmosphere. The scaling factor was chosen to be 
1/[1 + £?o/(87rP qs )], where P qs is the quiet-Sun pressure in the absence of the 
spot. This scaling factor works, although a more conservative scaling will be 
tested in the future. 

Above the model S atmosphere, the density and pressure fall rapidly and 
hydrostatic balance requires an extreme scaling of the Lorentz force. However 
since we are not aiming to realistically model waves that reach these heights 
(we only want to damp them) this is acceptable. Instead of requiring force 
balance in this purely artificial region, we have chosen to put Po(x,y,z) = 
P qs (z)Po(x,y,zo)/P cls (zo) where zq is the top of the model S atmosphere and 
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Pq S (z) is the pressure stratification of the system in the absence of the spot (the 
quiet-Sun value). The density was treated similarly. 

In this paper we follow Schluter and Temesvary (1958), Deinzer (1965), and 
Schussler and Rempel (2005) amongst others in concentrating on axisymmetric 
self-similar solutions. For this class of models, the vertical component of the 
magnetic field is assumed to satisfy 

B 0z (r, z)=B Q (ry/H(z)) H(z), (14) 

where the functions Q and H satisfy Q(0) = H(0) = 1 but are otherwise 
arbitrary functions, and Bo is a scalar measure of the vertical-field strength 
at the surface. The usual practice (Solanki, 2003), which we adopt in this paper, 
is to assume Q is a Gaussian, 

Q(r) = e-^ r2 / R o, (15) 

where i?o is the half-width at half-maximum. We chose i?o = 10 Mm to cor- 
respond to the observed value for sunspot AR9787, since Ro is the half- width 
at half-maximum of the model sunspot at the surface. The value of i? is fixed 
throughout this paper. For the function H we chose the exponential function 

H(z)=e z/a , (16) 

with a — 6.25 Mm. This choice is rather arbitrary and will certainly be varied 
in the near future. In this paper we concentrate on varying Bo, the peak mag- 
netic field at the surface. Having thus prescribed the formula for Bo z (r, z), we 
determine B 0r (r 7 z) by the requirement that V • B = 0. Once we have such a 
hydrostatic solution, we adjust the pressure everywhere to make it convectively 
stable in the manner described above. 

5. Comparison between Simulations and Observations 

In this section we want to compare the observed /-mode cross-covariance from 
AR9787 and the simulations. The computational domain is —40 Mm< x < 
105 Mm, -36.25 Mm< y < 36.25 Mm, and -12.5 Mm< z < 1.5 Mm. The 
sunspot axis is x = y = 0. The relationship of the computational box to the 
observations is shown by the white rectangle in Figure 4. 

5.1. The Initial Conditions for the Simulation 

In this section we discuss the choice of the initial conditions for the simula- 
tion. The aim is to allow a direct comparison between the simulated wave field 
(Section 4) and the observed cross-covariance function (Section 3). Since the 
observed cross-covariance uses the z component of the wave velocity as input 
data, we choose to use v' z from the simulations for the comparison (this is also 
in accordance with a deeper interpretation of the cross-covariance, see Campillo 
and Paul, 2003). The vertical component of velocity of any /-mode wave packet 
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Figure 6. The initial distribution of /-mode amplitudes used in the simulations, A^, as a 
function of kR@. 



propagating in the +x direction in a horizontally homogeneous atmosphere is of 
the form 

v' z (x,y,z,t) = ReJ2 A ke kz e iKx ~ xo) ~ iu ' kt ~ Jkt , (17) 
k 

where Ak are complex amplitudes, and exp(fcz) and u>k = \fgk~ are the /-mode 
eigenfunction and eigenfrcqucncy at wavcnumbcr k. We introduced the reference 
coordinate xq = —A, where A = 40 Mm is defined in the previous section as the 
distance between F and the sunspot. Given the evolution Equations (8) and (9), 
this wave packet is uniquely determined by the initial conditions 

v 1 = RcJ2(™ + z) A kC kz+ik{x - X0 \ (18) 

k 

£ = Rc^(-Z + \i)uj- 1 A k c kz+ik{x - X0 \ (19) 

k 

Thus our problem is reduced to fixing the amplitudes Ak- In this study, the 
amplitudes Ak are real (all waves are in phase at x — x n and t = 0) and are shown 
in Figure 6. This choice was simply the result of requiring that the simulated 
v z (x,y,z — 0,t) and the cross-covariance look approximately the same, in the 
far field and in the absence of the sunspot. As will be shown in the next section, 
this spectrum of wavcnumbers is sufficiently accurate for the present study. A 
more systematic analysis is planned. One important difficulty that arises is the 
existence of background solar noise which is not easy to model and has been 
ignored in Equation (17). 
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5.2. Magneto- Acoustic Waves 

First we consider a sunspot model with Bq = 3 kG. The simulation enables us to 
understand what happens below the solar surface. The theory of the interaction 
of solar waves with sunspots had been developed using ray theory and two- 
dimensional numerical simulations (e.g. Cally, 2007; Bogdan et al, 1996). We 
find what appears to be very similar physics in our three-dimensional simula- 
tions: the strong "absorption" of the /modes is consistent with partial conversion 
of the incoming waves into slow magnetoacoustic waves which propagate along 
the field lines. This can be seen in Figure 7. The downward propagating slow 
magnetoacoustic modes experience a decrease in their speed as they propagate 
downwards and are therefore shifted to increasingly short wavelengths. Mode 
conversion is a very robust feature of this type of simulation. We note that the 
slow magnetoacoustic waves are much easier to see in the x-component of wave 
velocity than in v' z , which is why Figure 7 shows the former. 

Since we are strongly damping short wavelengths, these magneto-acoustic 
modes rapidly decay. Any upward propagating wave encounters the damping 
buffer situated above the photosphere and is also damped. In principle this 
decay is unphysical in both instances, however since neither the downward nor 
upward propagating waves return to the surface this is not undesirable. 

5.3. Wave Forms 

In this section we compare the simulated v' z at z = with the observed cross- 
covariance. We emphasize that at the moment this comparison is not expected 
to be appropriate in the immediate vicinity of the sunspot. 

Figure 8 shows such a comparison between simulations and observations for 
times t — 40 minutes, t = 100 minutes and t = 130 minutes. The peak field 
strength used in this simulation was Bq — 3 kG and the radius Rq = 10 Mm. 
The comparison appears to be very good at time t = 130 minutes, at which time 
the wavepacket has completely traversed the sunspot. At this time, all aspects 
of the waveform would appear to have been approximately reproduced by the 
simulation: amplitudes (including "absorption"), phases, and spatial spectrum. 
The match appears to be less good for t = 40 minutes; this will be discussed 
below. 

To be more precise, however, the match between the simulations and obser- 
vations must be better quantified. In order to reduce the observational noise, we 
have averaged the cross-covariance function in the y-direction over two bands. 
Both bands, shown in Figure 9, are 14.5 Mm wide in the y-direction. The bands 
are labeled A and B and centered around y = and y = ±29 Mm. Band A 
is centered on the spot, while band B acts as a reference. The simulations are 
averaged in the same manner for comparison with the observations. 

In Figure 10 we have plotted both the simulation (v' z , thick red lines) and the 
observed cross-covariance (C, thin bluelines). In the top panels we compare the 
wave propagation at the edge of the computational box, i.e. in band B. This part 
of the wave is little affected by the sunspot and the match between the simulation 
and observation is quite good since the initial amplitude spectrum Ak was chosen 
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Figure 7. This plot shows pov' x in the x-z plane, through the sunspot axis. The upper frames 
are for time t = 34 min (before the /-mode wave packet crosses the sunspot) and the lower 
frames arc for t = 84 min (after). The black curves with equation B z (r,z) = B z (r = 0,z)/2 
give an estimate of the "width" of the sunspot. The conversion of the incoming / modes into 
slow magnctoacoustic modes is evident in the lower frames. 
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Figure 8. Comparison of the simulated vertical velocity and the observed cross-covariance. 
In each panel the upper frame shows the observed cross-covariance, the bottom panel the 
simulated wave packet. The panels correspond to t = 40 minutes, t = 100 minutes, and 
t = 130 minutes, from top to bottom. The black circles of radius Rq = 10 Mm indicate the 
location of the sunspot. 
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Figure 9. Sketch of the bands over which the data have been averaged in the y direction in 
later plots. Band A is centered on the spot, while band B acts as a reference. 



accordingly. In the lower panels of Figure 10 we see the results for the waves 
passing through the spot, i.e. band A. First we notice that the wave amplitude is 
remarkably well reproduced in the simulations, meaning that the observed wave 
"absorption" is consistent with mode conversion. For t = 130 minutes, after 
the waves have crossed the sunspot, there is no appreciable phase shift between 
the simulation and the observations. The match, however, is somewhat worse 
at t = 40 minutes and there is an obvious phase shift between the two wave 
forms with the observations lagging behind the simulations. Given the value of 
the phase shift, we suspect that it is due to the effect of the moat flow. The 
moat flow is a horizontal outflow from the sunspot, which we have measured 
by tracking the small moving magnetic features. The observed moat velocity 
(averaged over nine days) has a peak value of 230 ms _1 at a distance of 25 Mm 
from the center of the sunspot and vanishes at a distance of 45 Mm. The solar 
waves moving through the flow are first slowed down (against the flow) and 
later sped up again (with the flow). We have not modeled this effect yet. It is 
reassuring, however, to see that the Doppler shift caused by the moat appears 
to have disappeared at t = 130 minutes. 

Reiterating, it appears that, by using numerical simulations, seismic signa- 
tures can be followed through their passage across the spot. 
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Figure 10. The top panels show the simulation (v' z , thick red lines) and the observed cross- 
covariance (C, thin blue lines) averaged in the y-direction across band B (reference), at three 
different times (t). The bottom panels show the simulation (v' z , So = 3 kG, thick red lines) 
and observed cross-covariance (C, thin blue lines) averaged in the y-direction across band A, 
at three different times (t). The effect of the sunspot is very easily seen, by comparing with 
band B. 
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Figure 11. Simulated vertical velocity (thick red lines) and observed cross-covariance (thin 
bluelines) at t = 130 minutes, (a) For the band at the edge of the domain (large y). (b) For the 
band passing through the sunspot (y = 0). The three panels correspond to simulations using 
Bo = 2 kG, 2.5 kG, and 3 kG. It is seen that the simulations with Bo = 3 kG provide the best 
match in terms of both amplitude and phase. The stripes indicate the location of the sunspot. 
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6. Constraining B 

The solutions shown thus far have been for a model sunspot with a peak vertical 
field strength Bo = 3 kG at the surface z = 0. The match is good, which raises 
the question of whether other field strengths would match as well. The answer to 
this question can be seen in Figures 11 where we show the comparison between 
the simulations and observations for different field strengths, Bo = 2 kG and 
Bo = 2.5 kG. The comparison is made some time after the wave packet has 
passed through the sunspot, so the issue of the moat flow does not arise. At this 
stage we restrict ourselves to commenting that qualitatively the match is best, 
in phase and amplitude, for Bo — 3 kG; for the other values of Bo there is an 
apparent phase mismatch along y = and the amplitude of that part of the 
wave passing through the spot is not sufficiently damped. The quantification of 
the "goodness of fit" between the simulation and observations, which will allow a 
more accurate determination of the field strength, will be the subject of a future 
study. 



7. Discussion 

We have performed three-dimensional MHD simulations of waves propagating 
through sunspot models. The computations are set up in such a way as to allow 
comparing observed cross-covariances (except in the immediate vicinity of the 
sunspot). The parameters of the sunspot model can be chosen in such a way that 
its helioseismic signature is in good agreement with helioseismic observations of 
sunspot AR9787. 

A qualitative study using / modes has enabled us to place a constraint, Bq > 
3 kG, on the sunspot's near-surface field strength. The remaining differences 
reflect real differences between the model and the observed sunspot, such as the 
moat flow. 

The model atmosphere that we have constructed also supports p modes with a 
dispersion relation that is very close to that of the Sun. Obviously, the p modes 
in combination with the / modes should enable us to place more substantial 
constraints on the subsurface structure of the sunspot. In particular, it should be 
possible to simultaneously constrain the magnetic field strength Bo, the sunspot 
radius Ro, and the magnetic field inclination (controlled by the parameter a). 
Of course, the parametric sunspot model that we have considered in this paper 
is just one particular model: we plan to consider other types of sunspot models 
in the future. 

In summary, we believe that we have shown that the full-waveform modeling 
of sunspots is feasible. 

Acknowledgements SOHO is a project of international collaboration be- 
tween ESA and NASA. We are grateful to Manfred Schiissler for insightful 
discussions. 



cameron_rev.tex; 12/02/2008; 9:58; p. 20 



Hclioscismology of Sunspots 



References 

Birch, A., Kosovichev, A.G., Price, G.H., Schlottmann, R.B.: 2001, Astrophys. J. 561, L229. 
Bogdan, T.J., Cally, P.S.: 1997, Proc. Roy. Soc. London, Series A 453, 919. 
Bogdan, T.J., Hindman, B.W., Cally, P.S., Charbonncau, P.: 1996, Astrophys. J. 465, 406. 
Braun, D.C.: 1995, Astrophys. J. 451, 859. 

Braun, D.C., Duvall, T.L. Jr, LaBronte, B.J.: 1987, Astrophys. J. 319, L27. 

Cally, P.: 2000, Solar Phys. 192, 395. 

Cally, P.: 2007, Astron. Nachr. 328, 286. 

Cally, P., Bogdan, T.: 1993, Astrophys. J. 402, 721. 

Cally, P., Bogdan, T.: 1997, Astrophys. J. 486, L67. 

Cally, P., Goossens, M.: 2008, Solar Phys. in press. 

Cally, P., Bogdan, T., Zweibel, E.G.: 1994, Astrophys. J. 437, 505. 

Cally, P., Crouch, A., Braun, D.C: 2003, Mon. Not, Roy. Astron. Soc. 346, 381. 

Cameron, R., Gizon, L., Daiffallah, K.: 2007, Astron. Nachr. 328 313. 

Campillo, M., Paul, A.: 2003, Science 299, 547. 

Colin de Verdiere, Y.: 2006, http://fr.arxiv.org/abs/math-ph/0610043/. 

Christensen-Dalsgaard, J.,Dappcn, W., Ajukov, S.V., Anderson, E.R., Antia, H.M., Basu, S., 
Baturin, V.A., Berthomicu, G., Chaboycr, B., Chitre, S.M., Cox, A.N., Dcmarque, P., 
Donatowicz, J., Dzicmbowski, W.A., Gabriel, M., Gough, D.O., Gucnther, D.B., Guzik, 
J. A., Harvey, J.W., Hill, F., Houdek, G., Iglesias, C.A., Kosovichev, A.G., Leibachcr, J.W., 
Morel, P., Profhtt, C.R., Provost, J., Reiter, J., Rhodes, E.J. Jr, Rogers, F.J., Roxburgh, 
I.W., Thompson, M.J., Ulrich, R.K.: 1996, Science 272, 1286. 

Couvidat, S., Birch, A.C., Kosovichev, A.G.,: 2006 Astrophys. J. 640, 516. 

Crouch, A.D., Cally, P.S., Charbonneau, P., Braun, D.C, Desjardins, M.: 2005, Mon. Not. 
Roy. Astron. Soc. 363, 1188. 

Deinzer, W.: 1965, Astrophys. J. 141, 548. 

Duvall, T.L. Jr, Kosovichev, A.G., Murawski, K.: 1998, Astrophys. J. 505, L55 
Duvall, T.L. Jr, Jefferies, S.M., Harvey, J.W., Pomerantz, M.A.: 1993, Nature 362, 430. 
Gizon, L., Birch, A.C.: 2002, Astrophys. J. 571, 966. 

Goucdard, P., Stehly, L., Brenguier, F., Campillo, M., Colin de Verdiere, Y., Larose, E., 
Margerin, L., Roux, Ph., Sanchez-Sesma, F.J., Shapiro, N., Weaver, R.: 2008, Geophysical 
Prospecting, accepted. 

Hanasoge, S.M., Duvall, T.L. Jr: 2007, Astron. Nachr. 328 319. 

Hanasoge, S.M., Duvall, T.L. Jr, Couvidat, S.: 2007, Astrophys. J. 664 1243. 

Hanasoge, S.M., Larsen, R.M., Duvall, T.L. Jr, DeRosa, M.L., Hurlburt, N.E., Schou, J., Roth, 
M., Christensen-Dalsgaard, J., Lele, S.K.: 2006, Astrophys. J. 648, 1268. 

Hindman, B.W., Zweibel, E.G., Cally, P.: 1996, Astrophys. J. 459 760. 

Khomcnko, E., Collados, M.: 2006, Astrophys. J. 653, 739. 

Kosovichev, A.G., Duvall, T.L. Jr, Scherrer, P.H.: 2000, Solar Phys. 192, 159. 

Lindsey, C, Braun, D.C: 1997, Astrophys. J. 485, 895. 

Parker, E.: 1979, Astrophys. J. 230, 905. 

Rosenthal, C.S., Julien, K.A.: 2000, Astrophys. J. 532, 1230. 
Schiisslcr, M., Rempcl, M.: 2005, Astron. Astrophys. 441, 337. 

Schliiter, A., Temesvary, S.: 1958, in Lehnert, B. (ed.) Electromagnetic Phenomena in Cosmical 

Physics, IAU Symp. 6, Cambridge University Press, Cambridge, 263. 
Schunker, H., Cally, P.S.: 2006, Mon. Not. Roy. Astron. Soc. 372, 551. 
Solanki, S.: 2003, Astron. Astrophys. Rev. 11, 153. 

Spruit, H.C.: 1991, in Toomre, J., Gough, D.O., (eds), Challenges to Theories of the Structure 
of Moderate Mass Stars, Lecture Notes in Physics, Vol. 388, Springer- Verlag, Berlin, 121. 

Spruit, H.C., Bogdan, T.: 1992, Astrophys. J. 391, L109. 

Zhao, J., Kosovichev, A.G., Duvall, T.L. Jr: 2001, Astrophys. J. 557, 384. 

Zhao, J., Georgobiani, D., Kosovichev, A.G., Benson, D., Stein, R.F., Nordlund, A: 2007, 
Astrophys. J. 659, 848. 



cameron_rev.tex; 12/02/2008; 9:58; p. 21 



cameron_rev.tex; 12/02/2008; 9:58; p. 22 



